1.1 Εγκατάσταση βιβλιοθήκης

Εγκατάσταση του πρώτου και κύριου πακέτου (και ενός βοηθητικού)
(Αυτό γίνεται μια φορά σε κάθε Η/Υ - εάν έχει εγκατασταθεί δεν χρειάζεται να δοθεί ξανά η εντολή αυτή)

## CRAN repository:
install.packages("easypackages", repos = "http://cran.us.r-project.org")

## Github:
devtools::install_github("drsimonj/twidlr")

1.2 Φόρτωση βιβλιοθήκης

Φόρτωση του εν λόγω πακέτου

library(easypackages)

1.3 Εγκατάσταση βιβλιοθηκών

Με την εντολή packages('insert_package_name_here') γίνεται η εγκατάσταση των πακέτων που μας ενδιαφέρουν¹

Δοκιμάστε το

1.4 Φόρτωση βιβλιοθηκών

Εντολή για φόρτωση των πακέτων

libraries("car", "psych", "tidyverse", "broom", "twidlr", "leaps", "relaimpo", 
    "gvlma", "MASS")

1.5 Working Directory

Με την ακόλουθη εντολή, βλέπουμε σε ποιόν φάκελο δουλεύουμε:

getwd()

1.6 Αλλαγή του working directory

Ενώ με την εντολή αυτή, αλλάζουμε τον φάκελο στον οποίο δουλεύουμε και αποθηκεύουμε δεδομένα
[π.χ., setwd("E:/") εάν θέλουμε να δουλεύουμε στον σκληρό δίσκο Ε]

setwd("E:/Academic/Lessons/<c5><e9><e4><e9><ea><dc> <cc><e1><e8><de><ec><e1><f4><e1> <c2><ef><f4><e1><ed><e9><ea><de><f2>/Labs/Labs")

2.1 Εισαγωγή αρχείου

Τώρα ας εισάγουμε το αρχείο το οποίο περιέχει τα δεδομένα τα οποία μας ενδιαφέρουν να αναλύσουμε

Aegean <- readxl::read_excel("G:/Academic/Lessons/EMB/Labs/Labs/Aegean Islands.xlsx")

Κάθε αρχείο πρέπει να έχει ένα όνομα στην R
Ποιό είναι το όνομα του αρχείου μας;

2.2 Δομή των δεδομένων μας

Ας δούμε τα δεδομένα τα οποία εισάγαμε στην R

Aegean

και τι είδους δεδομένα είναι αυτά (ποιοτικά; ποσοτικά;)

str(Aegean)

## Classes 'tbl_df', 'tbl' and 'data.frame':    59 obs. of  14 variables:
##  $ Island   : chr  "Alonissos" "Amorgos" "Anafi" "Andros" ...
##  $ Area     : num  65 121 38 380 19.7 ...
##  $ Elevation: num  475 823 584 1003 378 ...
##  $ MAT      : num  111 133 139 130 139 ...
##  $ MAP      : num  489 498 498 393 362 ...
##  $ Temp     : num  161 175 181 170 181 ...
##  $ Rain     : num  495 528 523 463 402 ...
##  $ Dis      : num  4.76 10.33 21.14 9.29 1.7 ...
##  $ Dm       : num  41 117.8 207.1 55.6 112.2 ...
##  $ Geology  : num  5 5 8 6 3 1 5 4 1 10 ...
##  $ Total    : num  515 728 504 915 311 ...
##  $ Regional : num  21 75 45 58 28 8 24 52 34 96 ...
##  $ Endemics : num  18 59 35 47 25 6 15 36 17 30 ...
##  $ SIE      : num  0 1 0 3 1 0 0 1 0 1 ...

Μπορούσαμε να πάρουμε την πληροφορία αυτή με άλλον τρόπο;

2.3 Λιγότερα δεκαδικά

Με την ακόλουθη εντολή, θα εμφανίζονται πλέον λιγότερα δεκαδικά:

options(digits = 2)

3.1 Έλεγχος κανονικότητας

normality <- lapply(Aegean[, 2:14], shapiro.test)

normality

## $Area
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.3, p-value = 3e-15
## 
## 
## $Elevation
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.8, p-value = 2e-06
## 
## 
## $MAT
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.8, p-value = 1e-07
## 
## 
## $MAP
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.9, p-value = 2e-04
## 
## 
## $Temp
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.8, p-value = 3e-07
## 
## 
## $Rain
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.9, p-value = 5e-05
## 
## 
## $Dis
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.8, p-value = 5e-07
## 
## 
## $Dm
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.9, p-value = 0.01
## 
## 
## $Geology
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.9, p-value = 1e-04
## 
## 
## $Total
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.9, p-value = 3e-04
## 
## 
## $Regional
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.6, p-value = 2e-11
## 
## 
## $Endemics
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.5, p-value = 3e-13
## 
## 
## $SIE
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.2, p-value = 3e-16

3.2 Έλεγχος συσχέτισης

Ας ελέγξουμε τη συσχέτιση μεταξύ των μεταβλητών μας
Εάν κάποιες έχουν συσχέτιση \(\geq 0.8-0.9\), πρέπει να τις αφαιρέσουμε από τη μεταγενέστερη ανάλυση (Γιατί;)

corr.test(Aegean[, 2:14], method = "spearman")

## Call:corr.test(x = Aegean[, 2:14], method = "spearman")
## Correlation matrix 
##            Area Elevation   MAT   MAP  Temp  Rain   Dis    Dm Geology
## Area       1.00      0.75 -0.39  0.10 -0.42  0.07  0.15 -0.41    0.79
## Elevation  0.75      1.00 -0.33  0.07 -0.38  0.04  0.04 -0.32    0.67
## MAT       -0.39     -0.33  1.00  0.00  0.96  0.05 -0.23  0.46   -0.32
## MAP        0.10      0.07  0.00  1.00  0.05  0.99  0.40 -0.26    0.03
## Temp      -0.42     -0.38  0.96  0.05  1.00  0.09 -0.19  0.33   -0.36
## Rain       0.07      0.04  0.05  0.99  0.09  1.00  0.38 -0.22    0.01
## Dis        0.15      0.04 -0.23  0.40 -0.19  0.38  1.00 -0.27    0.22
## Dm        -0.41     -0.32  0.46 -0.26  0.33 -0.22 -0.27  1.00   -0.23
## Geology    0.79      0.67 -0.32  0.03 -0.36  0.01  0.22 -0.23    1.00
## Total      0.95      0.77 -0.42  0.09 -0.43  0.05  0.16 -0.47    0.79
## Regional   0.63      0.65 -0.19  0.11 -0.24  0.10  0.10 -0.17    0.58
## Endemics   0.43      0.45 -0.03 -0.15 -0.12 -0.15 -0.07  0.18    0.45
## SIE        0.62      0.53 -0.31  0.11 -0.32  0.09  0.35 -0.32    0.50
##           Total Regional Endemics   SIE
## Area       0.95     0.63     0.43  0.62
## Elevation  0.77     0.65     0.45  0.53
## MAT       -0.42    -0.19    -0.03 -0.31
## MAP        0.09     0.11    -0.15  0.11
## Temp      -0.43    -0.24    -0.12 -0.32
## Rain       0.05     0.10    -0.15  0.09
## Dis        0.16     0.10    -0.07  0.35
## Dm        -0.47    -0.17     0.18 -0.32
## Geology    0.79     0.58     0.45  0.50
## Total      1.00     0.67     0.48  0.67
## Regional   0.67     1.00     0.77  0.53
## Endemics   0.48     0.77     1.00  0.49
## SIE        0.67     0.53     0.49  1.00
## Sample Size 
## [1] 59
## Probability values (Entries above the diagonal are adjusted for multiple tests.) 
##           Area Elevation  MAT  MAP Temp Rain  Dis   Dm Geology Total
## Area      0.00      0.00 0.11 1.00 0.06 1.00 1.00 0.06    0.00  0.00
## Elevation 0.00      0.00 0.47 1.00 0.14 1.00 1.00 0.53    0.00  0.00
## MAT       0.00      0.01 0.00 1.00 0.00 1.00 1.00 0.02    0.54  0.06
## MAP       0.46      0.61 0.98 0.00 1.00 0.00 0.08 1.00    1.00  1.00
## Temp      0.00      0.00 0.00 0.68 0.00 1.00 1.00 0.47    0.25  0.04
## Rain      0.58      0.77 0.73 0.00 0.50 0.00 0.14 1.00    1.00  1.00
## Dis       0.27      0.75 0.08 0.00 0.15 0.00 0.00 1.00    1.00  1.00
## Dm        0.00      0.01 0.00 0.05 0.01 0.09 0.04 0.00    1.00  0.01
## Geology   0.00      0.00 0.01 0.84 0.01 0.95 0.10 0.08    0.00  0.00
## Total     0.00      0.00 0.00 0.50 0.00 0.69 0.22 0.00    0.00  0.00
## Regional  0.00      0.00 0.14 0.39 0.07 0.45 0.45 0.20    0.00  0.00
## Endemics  0.00      0.00 0.79 0.26 0.36 0.25 0.62 0.17    0.00  0.00
## SIE       0.00      0.00 0.02 0.39 0.01 0.50 0.01 0.01    0.00  0.00
##           Regional Endemics  SIE
## Area             0     0.04 0.00
## Elevation        0     0.02 0.00
## MAT              1     1.00 0.62
## MAP              1     1.00 1.00
## Temp             1     1.00 0.53
## Rain             1     1.00 1.00
## Dis              1     1.00 0.27
## Dm               1     1.00 0.53
## Geology          0     0.02 0.00
## Total            0     0.01 0.00
## Regional         0     0.00 0.00
## Endemics         0     0.00 0.01
## SIE              0     0.00 0.00
## 
##  To see confidence intervals of the correlations, print with the short=FALSE option

Διακρίνετε κάποιες μεταβλητές με τόσο μεγάλη συσχέτιση;

usdm::vifcor(Aegean %>% as.data.frame() %>% dplyr::select(Area:Geology), th = 0.8)

## 2 variables from the 9 input variables have collinearity problem: 
##  
## MAT MAP 
## 
## After excluding the collinear variables, the linear correlation coefficients ranges between: 
## min correlation ( Rain ~ Elevation ):  -0.021 
## max correlation ( Geology ~ Elevation ):  0.75 
## 
## ---------- VIFs of the remained variables -------- 
##   Variables VIF
## 1      Area 2.6
## 2 Elevation 3.9
## 3      Temp 1.9
## 4      Rain 1.3
## 5       Dis 1.7
## 6        Dm 1.2
## 7   Geology 2.4

3.3 Οπτικοποίηση αποτελεσμάτων

Μπορούμε να οπτικοποιήσουμε τα αποτελέσματα μας και να διαπιστώσουμε εάν κάποια μεταβλητή δεν εμφανίζει κανονική κατανομή

scatterplotMatrix(Aegean[, 2:14], spread = F, smoother.args = list(lty = 2), 
    main = "Scatter Plot Matrix")

3.4 Κανονικοποίηση των δεδομένων μας

Καθώς τα δεδομένα μας δεν έχουν κανονική κατανομή, θα πρέπει να τα κανονικοποιήσουμε. Αυτό σημαίνει ότι θα πρέπει να τα λογαριθμήσουμε. Συνεπώς, πρέπει να ελέγξουμε εάν κάποιες μεταβλητές περιέχουν το 0. Η ακόλουθη εντολή μας βοηθάει να το διαπιστώσουμε αυτό:

range(Aegean$SIE)

## [1]   0 172

Κάντε το και για τις υπόλοιπες μεταβλητές (Δηλαδή: range(Aegean$insert_IV_name_here))

3.5 Λογαρίθμηση

Εφ’όσον έχουμε ελέγξει όλες τις μεταβλητές μας, μπορούμε να τις λογαριθμήσουμε:

Aegean <- Aegean %>% mutate(Area = log10(Area), Total = log10(Total), Endemics = log10(Endemics), 
    Regional = log10(Regional), SIE = log10(SIE + 1))

Και εάν θέλουμε να λογαριθμήσουμε όλες τις μεταβλητές που μας ενδιαφέρουν, ένας γρήγορος τρόπος είναι ο ακόλουθος:

Aegean_log <- lapply(Aegean[, 2:13], as.numeric) %>% as_tibble() %>% log10() %>% 
    dplyr::mutate(Island = Aegean$Island, SIE = log10(Aegean$SIE + 1)) %>% dplyr::select(Island, 
    everything(), SIE)

3.6 Δομή των δεδομένων μας (ξανά)

Ας δούμε τα δεδομένα μας ξανά:

Aegean_log

4.1 Πολλαπλή γραμμική παλινδρόμηση με την μέθοδο enter

Τώρα είμαστε πλέον σε θέση να δημιουργήσουμε 4 μοντέλα πολλαπλής γραμμικής παλινδρόμησης, ένα για κάθε κατηγορία φυτικών taxa [Συνολικός αριθμός ειδών, Ενδημικά είδη, Περιφερειακά ενδημικά είδη και Τοπικά νησιωτικά είδη (Total, Endemics, Regional & SIE, αντίστοιχα)] σε σχέση με όλες τις ανεξάρτητες μεταβλητές μας². Θα χρησιμοποιήσουμε αρχικά την μέθοδο enter:

enter1 <- lm(Total ~ Area + Elevation + Rain + Temp + Dm + Dis + Geology, data = Aegean_log)

enter2 <- lm(Endemics ~ Area + Elevation + Rain + Temp + Dm + Dis + Geology, 
    data = Aegean_log)

enter3 <- lm(SIE ~ Area + Elevation + Rain + Temp + Dm + Dis + Geology, data = Aegean_log)

enter4 <- lm(Regional ~ Area + Elevation + Rain + Temp + Dm + Dis + Geology, 
    data = Aegean_log)

ΠΡΟΣΟΧΗ
Τα ανωτέρω μοντέλα δεν δείχνουν ποιές είναι οι μεταβλητές οι οποίες επηρεάζουν τα πρότυπα ποικιλότητας. Δείχνουν απλά τι ποσοστό εξηγείται εάν συμπεριλαμβάνονται όλες οι ανεξάρτητες μεταβλητές στο μοντέλο μας

Ας δούμε τώρα την περίληψη για κάθε ένα εξ αυτών των μοντέλων:

summary(enter1)

## 
## Call:
## stats::lm(formula = formula, data = data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.30188 -0.03583 -0.00143  0.04867  0.16530 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.73748    1.01870    1.71    0.094 .  
## Area         0.23040    0.03857    5.97  2.3e-07 ***
## Elevation    0.01921    0.07729    0.25    0.805    
## Rain         0.12354    0.11986    1.03    0.308    
## Temp         0.06476    0.41385    0.16    0.876    
## Dm          -0.03509    0.02071   -1.69    0.096 .  
## Dis          0.00981    0.02605    0.38    0.708    
## Geology      0.15085    0.06651    2.27    0.028 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.088 on 51 degrees of freedom
## Multiple R-squared:  0.854,  Adjusted R-squared:  0.834 
## F-statistic: 42.7 on 7 and 51 DF,  p-value: <2e-16

summary(enter2)

## 
## Call:
## stats::lm(formula = formula, data = data)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -0.501 -0.130  0.030  0.135  0.554 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)  -1.7112     3.0680   -0.56     0.58
## Area          0.1802     0.1161    1.55     0.13
## Elevation     0.3017     0.2328    1.30     0.20
## Rain         -0.3774     0.3610   -1.05     0.30
## Temp          1.2387     1.2464    0.99     0.33
## Dm            0.0821     0.0624    1.32     0.19
## Dis          -0.0205     0.0785   -0.26     0.79
## Geology       0.1331     0.2003    0.66     0.51
## 
## Residual standard error: 0.26 on 51 degrees of freedom
## Multiple R-squared:  0.383,  Adjusted R-squared:  0.298 
## F-statistic: 4.52 on 7 and 51 DF,  p-value: 0.000558

summary(enter3)

## 
## Call:
## stats::lm(formula = formula, data = data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.5765 -0.1688 -0.0096  0.1368  0.9108 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)  -0.4846     3.4535   -0.14   0.8890   
## Area          0.3987     0.1307    3.05   0.0036 **
## Elevation     0.4987     0.2620    1.90   0.0627 . 
## Rain         -0.1677     0.4063   -0.41   0.6816   
## Temp         -0.3976     1.4030   -0.28   0.7780   
## Dm           -0.0172     0.0702   -0.25   0.8070   
## Dis           0.1056     0.0883    1.20   0.2375   
## Geology      -0.1697     0.2255   -0.75   0.4553   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3 on 51 degrees of freedom
## Multiple R-squared:  0.602,  Adjusted R-squared:  0.547 
## F-statistic:   11 on 7 and 51 DF,  p-value: 2.18e-08

summary(enter4)

## 
## Call:
## stats::lm(formula = formula, data = data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.5901 -0.0956  0.0419  0.1268  0.4315 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept) -3.27429    2.48522   -1.32    0.194  
## Area         0.23149    0.09408    2.46    0.017 *
## Elevation    0.31825    0.18856    1.69    0.098 .
## Rain         0.22590    0.29241    0.77    0.443  
## Temp         1.25730    1.00964    1.25    0.219  
## Dm           0.02453    0.05053    0.49    0.629  
## Dis          0.00645    0.06356    0.10    0.920  
## Geology      0.10430    0.16226    0.64    0.523  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.21 on 51 degrees of freedom
## Multiple R-squared:  0.557,  Adjusted R-squared:  0.496 
## F-statistic: 9.15 on 7 and 51 DF,  p-value: 2.83e-07

4.2 Πολλαπλή γραμμική παλινδρόμηση με την μέθοδο της βηματικής πολλαπλής γραμμικής παλινδρόμησης

Με την μέθοδο αυτή θα μπορέσουμε να δούμε ποιές μεταβλητές εμπερικλείονται στο βέλτιστο μοντέλο για τον συνολικό αριθμό των ειδών

stepm1 <- stepAIC(lm(Total ~ Area + Elevation + Rain + Temp + Dm + Dis + Geology, 
    data = Aegean_log), scale = 0, direction = "both", trace = 1, keep = NULL, 
    steps = 1000, k = 2)

summary(stepm1)

Κάντε το ίδιο και για τις υπόλοιπες τρεις κατηγορίες ειδών (Endemics, SIE, Regional). Ποιές μεταβλητές εξηγούν και σε τι ποσοστό τα πρότυπα της ποικιλότητας για αυτές τις κατηγορίες;

ΠΡΟΣΟΧΗ: Μας ενδιαφέρει το adjusted R squared

4.3 Ποιό μοντέλο είναι καλύτερο;

Με την βοήθεια του Akaike Information Criterion (εν συντομία: AIC), μπορούμε να διαπιστώσουμε ποιό μοντέλο είναι “καλύτερο”. Όσο μικρότερη η τιμή του AIC, τόσο καλύτερο το μοντέλο.

AIC(enter1, stepm1)

Κάντε το ίδιο και για τις υπόλοιπες τρεις κατηγορίες ειδών (Endemics, SIE, Regional). Ποιό από τα δύο μοντέλα για κάθε κατηγορία ειδών είναι καλύτερο;

4.4 Δημιουργία μιας λειτουργίας για τον υπολογισμό της σχετικής σημαντικότητας των μεταβλητών

relweights <- function(fit, ...) {
    R <- cor(fit$model)
    nvar <- ncol(R)
    rxx <- R[2:nvar, 2:nvar]
    rxy <- R[2:nvar, 1]
    svd <- eigen(rxx)
    evec <- svd$vectors
    ev <- svd$values
    delta <- diag(sqrt(ev))
    lambda <- evec %*% delta %*% t(evec)
    lambdasq <- lambda^2
    beta <- solve(lambda) %*% rxy
    rsquare <- colSums(beta^2)
    rawwgt <- lambdasq %*% beta^2
    import <- (rawwgt/rsquare) * 100
    import <- as.data.frame(import)
    row.names(import) <- names(fit$model[2:nvar])
    names(import) <- "Weights"
    import <- import[order(import), 1, drop = FALSE]
    dotchart(import$Weights, labels = row.names(import), xlab = "% of R-Square", 
        pch = 19, main = "Relative Importance of Predictor Variables", sub = paste("Total R-Square=", 
            round(rsquare, digits = 3)), ...)
    return(import)
}

4.5 Υπολογισμός της σχετικής σημαντικότητας των μεταβλητών

Μετά την περισπωμένη (~) βάζουμε όποιες μεταβλητές έχουν βγει σημαντικές από την βηματική πολλαπλή γραμμική παλινδρόμηση:

rel.1 <- lm(Total ~ Area + Dm + Geology, data = Aegean_log)

relweights(rel.1, col = "blue")

Κάντε το ίδιο και για τις υπόλοιπες τρεις κατηγορίες ειδών (Endemics, SIE, Regional). Ποιά μεταβλητή είναι σημαντικότερη για την εκάστοτε κατηγορία ειδών;

4.6 Πολλαπλή γραμμική παλινδρόμηση με την μέθοδο των ολικών υποσυνόλων

Μέχρι στιγμής, έχουμε δει δύο διαφορετικές μεθόδους, όσον αφορά την πολλαπλή γραμμική παλινδρόμηση. Όμως, δεν μπορούμε να είμαστε 100% σίγουροι ότι έχουμε φθάσει στο ορθό συμπέρασμα, καθότι ένα από τα μειονεκτήματα της βηματικής πολλαπλής παλινδρόμησης, είναι ότι μπορεί να “κολλήσει” σε ένα υποβέλτιστο πλατώ. Τον σκόπελο αυτό, μπορούμε να τον αποφύγουμε χρησιμοποιώντας την μέθοδο των ολικών υποσυνόλων:

bes.tot <- regsubsets(Total ~ Area + Elevation + Rain + Temp + Dm + Dis + Geology, 
    data = Aegean_log, nbest = 4)

Κάντε το ίδιο και για τις υπόλοιπες τρεις κατηγορίες ειδών (Endemics, SIE, Regional).

4.7 Απεικόνιση των αποτελεσμάτων του προηγούμενου βηματος

Το γράφημα αυτό θα μας δείξει ποιός είναι ο πραγματικά βέλτιστος συνδυασμός των ανεξάρτητων μεταβλητών:

plot(bes.tot, scale = "adjr2")

Κάντε το ίδιο και για τις υπόλοιπες τρεις κατηγορίες ειδών (Endemics, SIE, Regional).

tot <- lm(Total ~ Area + Rain + Dm + Geology, data = Aegean_log)

summary(tot)

## 
## Call:
## stats::lm(formula = formula, data = data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.30141 -0.03724  0.00481  0.04501  0.16857 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   1.9077     0.3196    5.97  1.9e-07 ***
## Area          0.2345     0.0301    7.79  2.1e-10 ***
## Rain          0.1341     0.1126    1.19    0.239    
## Dm           -0.0364     0.0197   -1.85    0.069 .  
## Geology       0.1543     0.0636    2.43    0.019 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.085 on 54 degrees of freedom
## Multiple R-squared:  0.854,  Adjusted R-squared:  0.843 
## F-statistic: 78.9 on 4 and 54 DF,  p-value: <2e-16

relweights(tot, col = "blue")

Κάντε το ίδιο και για τις υπόλοιπες τρεις κατηγορίες ειδών (Endemics, SIE, Regional). Ποιές είναι οι σημαντικότερες μεταβλητές;

4.8 Έλεγχος των αποτελεσμάτων της παλινδρόμησης

Τώρα θα ελέγξουμε κατ’αρχάς εάν οι μεταβλητές του βέλτιστου μοντέλου μας παραβιάζουν μια από τις αρχές της πολλαπλής γραμμικής παλινδρόμησης:

sqrt(vif(tot)) > 2

##    Area    Rain      Dm Geology 
##   FALSE   FALSE   FALSE   FALSE

Κάντε το ίδιο και για τις υπόλοιπες τρεις κατηγορίες ειδών (Endemics, SIE, Regional).

4.9 Έλεγχος πιθανοτήτων

Στο σημείο αυτό θα χρειαστεί να ελέγξουμε εάν η πιθανότητα κάθε μεταβλητής η οποία εμπερικλείεται στο τελικό μας μοντέλο είναι πραγματικά στατιστικώς σημαντική:

##--------------------------------------------------------------------------
## Let's have a look on the p-values for each of the IVs that have entered
## our final model
##--------------------------------------------------------------------------
summary(tot)

## 
## Call:
## stats::lm(formula = formula, data = data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.30141 -0.03724  0.00481  0.04501  0.16857 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   1.9077     0.3196    5.97  1.9e-07 ***
## Area          0.2345     0.0301    7.79  2.1e-10 ***
## Rain          0.1341     0.1126    1.19    0.239    
## Dm           -0.0364     0.0197   -1.85    0.069 .  
## Geology       0.1543     0.0636    2.43    0.019 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.085 on 54 degrees of freedom
## Multiple R-squared:  0.854,  Adjusted R-squared:  0.843 
## F-statistic: 78.9 on 4 and 54 DF,  p-value: <2e-16

##--------------------------------------------------------------------------

##--------------------------------------------------------------------------
## Inside the parenthesis we include the above p-values
##--------------------------------------------------------------------------
p <- c(0, 0.239, 0.069, 0.019)
##--------------------------------------------------------------------------


##--------------------------------------------------------------------------
## Level of significance
##--------------------------------------------------------------------------
alp <- 0.05
##--------------------------------------------------------------------------


##--------------------------------------------------------------------------
## If the model's p-values are equal or smaller than the p-values we obtain
## from the following command, then they are ok
##--------------------------------------------------------------------------
p.adjust(p, method = "bonferroni", n = length(p))

## [1] 0.000 0.956 0.276 0.076

##--------------------------------------------------------------------------

Κάντε το ίδιο και για τις υπόλοιπες τρεις κατηγορίες ειδών (Endemics, SIE, Regional).